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Abstract. We present the new open-source spherically-symmetric general- 
relativistic (GR) hydrodynamics code GRID. It is based on the Eulcrian formulation 
of GR hydrodynamics (GRHD) put forth by Romero-Ibancz-Gourgoulhon and 
employs radial-gauge, polar-slicing coordinates in which the 3+1 equations 
simplify substantially. We discretize the GRHD equations with a finite- 
volume scheme, employing piecewise-parabolic reconstruction and an approximate 
Riemann solver. GRID is intended for the simulation of stellar collapse to neutron 
stars and black holes and will also serve as a testbed for modeling technology to 
be incorporated in multi-D GR codes. Its GRHD part is coupled to various finite- 
temperature microphysical equations of state in tabulated form that we make 
available with GRID. An approximate deleptonization scheme for the collapse phase 
and a neutrino-leakagc/hcating scheme for the postbounce epoch are included and 
described. We also derive the equations for effective rotation in ID and implement 
them in GRID. We present an array of standard test calculations and also show 
how simple analytic equations of state in combination with presupernova models 
from stellar evolutionary calculations can be used to study qualitative aspects of 
black hole formation in failing rotating core-collapse supernovae. In addition, we 
present a simulation with microphysical EOS and neutrino leakage/heating of a 
failing core-collapse supernova and black hole formation in a presupernova model 
of a 40- Mq zero-age main-sequence star. We find good agreement on the time of 
black hole formation (within 20%) and last stable protoneutron star mass (within 
10%) with predictions from simulations with full Boltzmann neutrino radiation 
hydrodynamics. 



PACS numbers: 04.25.D-,04.40.Dg,97.10.Kc,97.60.Bw,97.60. Jd,97.60.Lf,26.60.Kp 
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1. Introduction 



Stellar core collapse is among the most energetic phenomena in the modern universe 
and liberates of the order of a few hundred [B]ethe (1 B = 10 51 erg) of gravitational 
energy as the core of a massive star (zero-age main-sequence [ZAMS] mass 8 — 10 Mq < 
M < 100 Mq) is compressed from a radius of ~ 1500 km and central density 
p c ~ 10 10 gcm -3 to ~ 15 km and p c in excess of nuclear density. Most (~ 99%) 
of this energy is ultimately radiated in neutrinos, but a small fraction (~ 1 B) may 
be converted into kinetic and internal energy of an outgoing shock wave and may 
result in a corc-collapse supernova explosion within the first seconds after collapse. 
The precise mode of conversion, the core-collapse supernova mechanism, is uncertain 
and has been the enigma of supernova theory for the past five decades (e.g., 0- 
0|). At the densities and velocities encountered in stellar collapse, the inclusion of 
general rclativistic effects is not an optional model sophistication, but a necessity for 
quantitatively and qualitatively reliable results. Importantly, general relativity (GR) 
predicts that the protoneutron star (PNS) formed in the initial collapse will undergo 
a second gravitational instability and collapse to a black hole (BH), if continued 
accretion pushes it over the maximum mass supported by the strong force and nucleon 
degeneracy. This may happen if the supernova mechanism fails and no explosion 
is launched or due to fallback accretion if an explosion occurs, but is too weak to 
unbind the entire stellar envelope @. In both cases, and provided sufficient angular 
momentum and its appropriate distribution in the progenitor star, the newly formed 
collapsar may become the central engine for a long-soft gamma-ray burst (GRB) 

II Ei. 



General relativistic computational models of stellar collapse have a long pedigree, 
starting with the spherically-symmetric (ID) Lagrangian work of May & White 
in the mid-1960s jl2j . based on the comoying GR hydrodynamics formulation in 
orthogonal coordinates by Misner & Sharp [13J and using a finite-difference scheme 
with an artificial viscosity [IH approach to handle shocks. Much subsequent ID 
GR work fl5T]2o| was based on this or similar approaches, including full radiation- 
hydrodynamics stellar collapse and core-collapse supernova simulations with finite- 
temperature microphysical equations of state (EOS) [21 - 25 1 . Eulerian formulations, 
more suited for extension to multi-D simulations, were introduced later and used 
maximal slicing f26T- l29| |. or radial-gauge, polar-slicing (RGPS) [1^. These schemes, 
with the exception of [30(, who employed pseudospectral methods, still used artificial 
viscosity approaches to shock treatment. More accurate, high-resolution shock- 
capturing (HRSC) approaches to GR stellar collapse based on higher-order Gudonov 



schemes and Ricmann solvers were introduced by Marti et al. [31[ and Yamada [32| 
in the Lagrangian context, by Marti et al. |33| in the fixed-background Eulerian case, 
and by Romero et al. and Noble [35| in the RGPS Eulerian frame. Yamada's 
approach was later extended to include microphysical EOS and radiation transport 
(3a . [37j . Gourgoulhon & Haensel (38j included an approximate neutrino transport 
treatment in their code. Preliminary results of Romero's code with a microphysical 
EOS and a neutrino leakage scheme were published in [3^, H(J ■ 

State-of-the-art simulations of stellar collapse and of the postbounce supernova 
evolution strongly suggest that multi-D dynamics is crucial for the core-collapse 
supernova mechanism to succeed in massive stars (e.g., I2|, 4, 4lM44|). Present multi-D 
core-collapse supernova codes are either Newtonian [44j-|46j or employ Newtonian 
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dynamics with relativistic corrections to the gravitational potential [!, 0, Q • Multi-D 
simulations in conformally-flat 47| or full GR traditionally relied on simple analytic 
EOS and polytropic initial models and neglected crucial neutrino effects (see, e.g., 



48W5l|). Only recently have the first axisymmetric (2D) [52|, |53j and 3D [5 



GR core collapse simulations become available that employ microphysical EOS and 
an approximate treatment of dclcptonization in the collapse phase, but postbounce 
neutrino transport, cooling, and heating are still not taken into account in these 
models. However, very recently, Miiller [5ri| has succeeded in implementing the 
complex and computationally-intensive radiation-transport scheme of 57[ in the 2D 
conformally-flat GR framework of [3 [4i| and first results are forthcoming 58 1 . 

In this article, we lay the foundations for a new and open approach to the stellar 
collapse and core-collapse supernova problem in GR. We discuss the formulation and 
implementation of the code GRID, a new, spherically-symmetric Eulcrian GR code for 
stellar collapse to neutron stars and black holes with approximate pre- and postbounce 
neutrino treatment. We release GRID and all its microphysics and input physics as 
open source to be downloaded from jhttp : //www. stellarcollapse . org. It is meant 
to complement open-source 3D GR codes such as Whisky [H^] that do not come with 
microphysics and neutrino approximations. At the same time, we intend GRID to 
serve as an efficient ID GR testbed for new modeling technology to be eventually 
incorporated in multi-D codes. In addition, GRID and its microphysics components can 
readily be adapted for use in the computational modeling of problems involving some 
or much of the same physics as in the stellar collapse problem, e.g., the postmerger 
phase of double neutron-star or black-hole - neutron-star coalescence. 

We base GRID on the conceptually simple and computationally efficient RGPS 
formalism of (30| . GRID, like the code of [jyjl, employs a Eulerian formulation of 
GR hydrodynamics with HRSC and works on non-equidistant grids. For the first 
time in the ID GR context, we derive and implement in GRID an extension of the 
ID GR hydrodynamics equations to include rotation in an effective fashion. For 
completeness and comparison of Newtonian and GR dynamics, GRID also implements 
ID Newtonian hydrodynamics. GRID operates with analytic EOS as well as with 
tabulated microphysical EOS through a general EOS interface. We discuss and provide 
EOS tables for the EOS of Lattimer-Swesty [(30| and the one of H. Shen et al. jpjs3|. 
Furthermore, we discuss and include in GRID the deleptonization treatment of [63j for 
the collapse phase and a postbounce 3-flavor neutrino treatment based on the leakage 
schemes of [64l . [65J as well as an approximate way of including neutrino heating. 

Due to these approximations in the neutrino treatment, GRID in its present form 
cannot be used for accurate simulations addressing the core-collapse supernova 
mechanism or neutrino-induced nucleosynthesis. However, we find that with the 
present treatment, GRID reproduces very well qualitatively the salient features of the 
postbounce evolution of core-collapse supernovae as predicted by full ID radiation- 
hydrodynamics simulations. Moreover, we find that GRID may be used to make 
quantitatively reliable predictions on the time of black hole formation in failing core- 
collapse supernovae and on the maximum mass of the PNS. 

This article is structured as follows. In section [3J we discuss our ID GR 
hydrodynamics and curvature equations and their implementation in GRID. Section [3] 
introduces the EOS provided with GRID and in section [4] we detail our pre-bounce 
deleptonization and postbounce leakage and neutrino heating schemes. A number of 
code tests and example simulations are presented in section [5] and |6] We wrap up and 
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conclude in section [7J 

We assume spacclikc signature ( — ,+,+,+) and, unless mentioned otherwise, use 
units of G = c = Mq = 1, but use cgs units for the microphysics and neutrino 
leakage/heating quantities. 



2. ID GR Hydrodynamics and Curvature Equations 



2.1. Curvature Equations in ID RGPS 

We follow 3(| 34 1 who formulate the 3 + 1 GR curvature and hydrodynamics equations 
in RGPS coordinates. In these coordinates and in spherical symmetry, the shift vector 
vanishes and the metric is diagonal and closely resembles the Schwarzschild metric. 
The invariant line element is 

ds 2 = g iiv dx p 'dx v , 

= - a(r, tfdt 2 + X(r, tfdr 2 + r 2 dVL 2 , (1) 

where a and X can be written more conveniently as functions of a metric potential, 
t), and the enclosed gravitational mass M grav (r, t) = m(r,t), 

-1/2 



a(r, t) = cxp [$(r, t)] 



X(r,t) 



1 



2m(r, t) 



(2) 



We assume ideal hydrodynamics for which the fluid stress-energy tensor and the matter 
current density are- 



= phu^u" + Pg^ v and J' 1 = pu* 



(3) 



where p is the baryonic density, P is the fluid pressure, h is the specific enthalpy equal 
to 1 + £ + P I p with e being the specific internal energy, is the four- velocity and, in 

— 1/2 

ID without rotation, is equal to [W/a, Wv r , 0, 0]. W = (l — w 2 ) is the Lorentz 
factor and v = Xv r . The equation for the gravitational mass needed for determining 
the metric coefficient X{r, t) of @ is derived from the Hamiltonian constraint equation 
and reads 



m(r, t) = 4tt [ (phW 2 -P + T^)r' 2 dr' 
Jo 



(4) 



Here, T^ n is the contribution to the gravitational mass from the energy and pressure 
of trapped neutrinos (see section FO)) . The expression for the metric potential $(?', t) 
is determined via the momentum constraints, taking into account the polar slicing 
condition that imposes tv K = K^, where K^j is the extrinsic curvature tensor (see 
30l . 35 1 for details). It reads, 

m(r', t) 



*(r,t) 



X s 



4irr'(phW 2 v 2 + P + r| 



dr' + $n 



(5) 



where analogous to (jlj), accounts for the effect of trapped neutrinos. $c 



is determined by matching the solution at the star's surface (r 
Schwarzschild metric, 



$(iZ*,t) = hx[a{R in t)] = iln 



1 - 



2m{R i ,,t) 



R+) to the 



(6) 



We use standard 2 nd order methods to perform the integrals in @ and (O and obtain 
values at cell centers as well as at cell interfaces. 
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2.2. GR Hydrodynamics m ID RGPS 



The evolution equations for the matter fields are derived from the local conservation 
laws for the stress-energy tensor, V ^T^" = 0, and for the matter current density 
V^J' 1 = 0. We write the GR hydrodynamics equations along the lines of the flux- 
conservative Valencia formulation (e.g., [66^68|) with modifications for spherically- 
symmetric flows proposed by [34| and neutrino sources. Derivation details are 
presented in |Appcndix A| 
We write the set of evolution equations as, 



X 



-F 



S, 



(7) 



where U is the set of conserved variables, F is their flux vector, and S is the vector 
containing gravitational, geometric, and neutrino-matter interaction sources and sinks. 
In ID and without rotation, U = [D, DY e , S r , r]. The conserved variables are functions 
of the primitive variables p, Y e , e, v, and P and are given by 

D = aXJ 1 = XpW , 
DY e = aXY e J f = XpWY e , 
S r = aXT tr = phW 2 v , 
r = a 2 T u - D = phW 2 



P - D 



(8) 

where Y e is the electron fraction, the number of electrons per baryon, and the only 
compositional variable needed to describe matter in nuclear statistical equilibrium 
(NSE). Note that there is a misprint in the central part of Eq. 9 of [34[ which 
is missing a factor of X which we have corrected here. The flux F is given by 
F = [Dv, DY e v, S r v + P,S r — Dv] and the sources and sinks are given by 



S 



2aP 



t - D)aX ( 8irrP 



aPX- 



Xr 



-iiaM 



(9) 



The source and sink terms Ry , Q'gV , Q V gr ,Qt' E i an d Qr' M are associated with 
neutrinos and are discussed in section U and derived in |Appcndix A| 

We use a semi-discrete approach and first discretize ([7]) in space, then apply the 
method of lines (MoL, [(3j|) and perform the time integration of the conserved variables 
via standard 2 nd or 3 rd order Runge-Kutta integrators with a Courant factor of 0.5. 



The spatial discretization follows the finite-volume approach (e.g., [3J, |68|) and all 
variables are defined at cell centers i and must be reconstructed (i.e., interpolated) 
at cell interfaces, where inter-cell fluxes are computed. This interpolation must be 
monotonic to ensure stability. We use the nominally 3 rd order (in smooth parts 
of the flow) piecewise-parabolic method (PPM, [7(3]) to interpolate the primitive 
variables and then set up the conserved variables at the cell interfaces. We 
also implement piecewise-constant reconstruction as well as piecewise-linear (total- 
variation-diminishing [TVD]) reconstruction with Van Leer's limiter 71(. The latter 
we use exclusively in the innermost 3 to 5 zones to avoid oscillations near the origin. 

Once the variables have been reconstructed at the cell interfaces, we evaluate the 
physical interface fluxes F i+1 / 2 with the HLLE Riemann solver (72|. The right-hand- 
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side (RHS) flux update term for Ui then reads. 



RHSi 



1 



r 2 Ar, 



X,. 



F 



i+l/2 



^-1/2^-1/2 
X~. 



1/2 



(10) 



-1/2 -A-i-1/2 

Gravitational, geometrical, and neutrino- matter interaction sources/sinks are not 
taken into account in the flux computation and are coupled into the MoL integration. 

After the update of the conserved variables D, DY e , S r and t, primitive variables p, 
Y e , v, e, and P{p, e, Y e ) must be extracted since they are needed for the next timestep. 
In the general case, the primitive variables (with the exception of Y e ) cannot be 
expressed algebraically in terms of the conserved variables (see, e.g., [67|). Hence, we 
employ an iterative approach and make an initial guess using P Q \d from the previous 
timestep, 

S r D t + D + P old (l - W 2 ) 
V = t + D + P m > P= XW> 6 = pW 11 

where we note that X can be calculated from the conserved variables as phW 2 — P = 
t + D. W is calculated from the estimate of v. We then call the EOS to obtain a new 
pressure and iterate this process using a Newton-Raphson method until convergence 
(we typically stop the iteration at a fractional pressure difference of 10~ 10 between 
iteration steps). 



2.3. Extension to 1.5D: Including Rotation 



Lagrangian spherically-symmetric stellar evolution codes have long included rotation 
and rotational effects in an approximate fashion (e.g., 73-75j|). The way this is 
typically done is to make the assumption that the star has constant angular velocity on 
spherical shells. In order to compute the effective specific centrifugal force acting on a 
fluid parcel, we compute the angular average of {Co x r) 2 on a spherical shell of radius r, 
which leads to /cent = 2/3 uj 2 r. In Newtonian Lagrangian calculations, specific angular 
momentum j = tor 2 is conserved by construction and the effective centrifugal force 
appears in the momentum equation. Relatively recently, such a n ap proach has also 
been taken in the Newtonian ID core collapse calculations of 7^, ZZl in order to take 



into account the effect of rotation approximately In the Eulcrian frame and in GR 
the situation is more complicated. We must solve an equation for angular momentum 
conservation on top of taking into account a centrifugal force term in the momentum 
equation. We begin by defining an azimuthal Eulerian velocity v^(= uS) and, in 



order to obtain a quantity of dimension velocity, we also define tv 



(note that 



= Wv v /r). With finite u*, T r * is finite and W become s W = {1-v 2 - 2/'iv 2 p )- 1 l 2 
in our effective approach. We provide derivation details in |Appendix A.2| and present 
here only the results. The modified stress-energy tensor leads to an additional equation 
for angular momentum conservation analogous to ([7]). 

1 / ar 2 \ 

dt(S <p ) + -d r (—F 4> )=S^ (12) 



where 



Stf, = phW 2 v v r , 

F^ = phW 2 v v rv = S^v 



= phW 2 avv lp X 



Anr z P 



(13) 
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Also, an additional term, accounting for the centrifugal force, 

appears on the RHS of the equation for S r . Finally, the change of the stress-energy 
tensor also has an effect on the metric potential $, whose equation is now given by 



— + Airrl phW\v z + -<) + P + t% 



(15) 



We implement this 1.5D treatment of rotation in GRID, but keep the metric diagonal. 
The 1.5D treatment should be rather accurate for slow rotation, and, as shown by fl7\ , 
will still capture qualitatively the effect of centrifugal support due to rapid rotation. 
For completeness, we note that the total angular momentum of the system (see, e.g., 
[78j |) is given by, 

J = Tl^gd 3 x = — / phXW 2 rv v r 2 dr , (16) 
Jo 3 Jo 

where we include a factor of 2/3 to account for the angular average. The rotation 
parameter /?, defined as the ratio T/\ W gra v| of rotational kinetic to gravitational energy 
is 

T/\W^ V \ = — — ^ — - , (17) 

where 

1 r°° Air c°° 

T=2j Q uT^y/=gePx = YJ Q P hXw2v l r2dr . ( 18 ) 

where again a factor of 2/3 in the last step is from performing an angular average. 

-^proper 

is given by, 

/•DO 

M propcr = 4tt / (p + pe) XWr 2 dr , (19) 
Jo 

and M grav is specified by (j4]). 



3. Equations of State (EOS) 



An EOS is needed to close the system of GR hydrodynamics equations and provide 
the pressure as well as other thermodynamic quantities as a function of density, 
temperature (or specific internal energy), and composition. In GRID, we include for 
test simulations the standard analytic polytropic (isentropic "cold", P = Kp T ) and 
the T-law EOS ("hot", P = (T — l)/?e). These are inappropriate for stellar collapse 
since they do not capture the stiffening of the EOS at nuclear density. An analytic 
EOS, able to cap ture this effect qualitatively and include nonisentropic effects, is the 
hybrid EOS [79( which we include in GRID and discuss in section 13.11 For a more 
realistic description of the thermodynamics of nuclear matter, an EOS built from a 
microphysical finite-temperature model for nuclear matter is needed. This is also a 
prerequisite for any kind of neutrino treatment, since crucial compositional information 
as well as chemical potentials must be derived from a microphysical model. Such 
microphysical EOS arc too complicated to be computed on the fly in a simulation and 
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are used in tabulated form with interpolation. GRID is able to handle such EOS and 
we provide tables at http://www.stellarcollapse.org/microphysics for the EOS 



of Lattimer & Swesty (|60(, LS EOS) and for the one of H. Shcn et al. ([6l|,l62j, HShcn 
EOS). The details of these tables and the routines facilitating their use are discussed 
in sections 13.21 and 13.31 



3.1. Hybrid EOS 

The hybrid EOS found widespread use in early multi-D simulations of rotating core 
collapse (e.g., (HHH), but was shown by [52], |53j to lead in some cases to qualitatively 



incorrect results for the collapse dynamics and the resulting gravitational wave signal. 
We include it in GRID, because its analytic nature provides for very fast calculations, 
allowing us to readily test the GR hydrodynamics of GRID. 
The hybrid EOS splits the pressure into a polytropic (cold) and a thermal component, 

P = Pcold + ^thermal . (20) 

The cold part is piecewise polytropic. It is composed of a polytropic EOS with 
r = Y\ for densities below nuclear (p nuc ) and another polytropic EOS with L = L 2 
for densities above p nuc . The two arc smoothly matched at p nuc which makes the 
polytropic constant K2 of the high-density part a function of the two Ts, of Ki, and 
of the transition density p nuc (see, e.g. |79j-|81| for a description of the procedure 
and detailed expressions). The thermal part is modeled via a L-law with F t h- It 
becomes relevant only after core bounce when shocks are present, making the flow 
nonadiabatic. Its contribution is determined via the thermal specific internal energy 
which is the difference between the primitive variable e and the cold specific internal 
energy, e t h = e - e co id- 
For collapse simulations, we set K x = 1.2435x 10 15 (Y e ) 4 ^ 3 [cgs] (the value appropriate 



for a relativistic degenerate gas of electrons, |80|,|82j) with Y e = 0.5. We choose a value 
below, but close to 4/3 for Li and typically set L 2 = 2.5 to mimic the stiff nuclear EOS 
above p nuc which we set to 2 x 10 14 gcm -3 . L t h we normally keep at 1.5 to model a 
mixture of relativistic (r = 4/3) and nonrclativistic (r = 5/3) thermal contributions. 
This leads to rapid shock propagation and explosion. When simulating BH formation 
with the hybrid EOS, we set I\h to smaller values. This reduces the postshock thermal 
pressure and leads to shock stagnation. 

3.2. Lattimer-Swesty EOS 

The LS EOS [(3(| is derived from a finite-temperature compressible liquid-droplet 
model [83d with a Skyrme nuclear force, uses the single heavy nucleus approximation, 
and assumes nuclear statistical equilibrium (NSE). NSE holds at T > 0.5 MeV which 
in core collapse and supernova matter is typically the case at p > few x 10 7 gcm~ 3 . 

The LS EOS routines are open source and available from the Stony Brook grourJJ. 
We employ their baryonic parts to generate tables with nuclear incomprcssibilitics 
K of 180 MeV, 220 MeV, and 375 MeV (the larger K , the stiffer the nuclear EOS). 
Hereafter, we refer to these Kq- variants of the LS EOS as LS180, LS220, and LS375. 
The symmetry energy S v is set in all variants to 29.3 MeV for all Kp. Electrons and 
photons are added using the routines provided by Timmes's EOS^ (84|. 



http: //www. astro . sunysb . edu/dswesty/lseos .html 
http: //cococubed . asu . edu/code_pages/eos . shtml 
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We compute the maximum cold neutron star masses for the three LS EOS variants 
by setting T = 0.1 MeV and assuming neutrino-less /3-equilibrium. The results 
are 1.83 M Q (2.13 M ), 2.04 M Q (2.41 -Mq), 2.72 M (3.35 M Q ) for gravitational 
(baryonic) mass and for K = 180 MeV, K = 220 MeV, and K a = 375 MeV, 
respectively. The coordinate radii of these maximum mass stars are 10.1 km, 10.6 km, 
and 12.3 km. 

Our LS EOS tables have 18 evenly-spaced points per decade in log 10 p ranging 
from 10 3 — 10 16 gcm -3 , 30 points per decade in log 10 T ranging from 10~ 2 — 
10 2 ' 4 MeV, and 50 points equally spaced in electron fraction from 0.035 to 0.53. 
This table resolution is sufficiently good to allow the use of simple and fast tri-linear 
interpolation (in log 10 (p), log 10 (T), Y e ), in collapse simulations while maintaining 
good thermodynamic consistency. In tests of adiabatic collapse, the inner-core entropy 
is conserved to ~ 1% from the onset of collapse to core bounce. 

To generate the LS EOS tables, we employ the LS EOS at densities above 10 8 gem" 3 , 
but, due to unreliable convergence, use linear extrapolation of the Helmholtz free 
energy F in Y e for Y e > 0.5 and in T at T < 0.06 MeV. Note that the latter is far away 
from NSE, but is never reached by core collapse tra ject ories at p > 10 8 gcm -3 . At 
densities below 10 8 gcm -3 , we use the Timmes EOS [84( and assume that the matter 
is an ideal gas composed of electrons, photons, neutrons, protons, alpha particles, and 
heavy nuclei with the average A and Z given by the LS EOS at the transition. 

Since the specific internal energies returned by the baryonic part of the Timmes 
EOS do not contain the nuclear binding energy, we shift the zero point of the Timmes 
EOS so that the returned specific internal energies are consistent with the LS EOS 
values at the transition point. For simplicity, we keep baryonic compositional variables 
fixed at the values obtained from the LS EOS at the transition density. These 
particular choices for the baryonic component have little effect at low densities where 
the thermodynamics are dominated by electrons at low to intermediate temperatures 
and by photons at high temperatures. However, for full core-collapse supernova 
simulations that intend to address also nuclear burning and nucleosynthesis aspects, 
a more involved consistent NSE/non-NSE EOS treatment involving the advection of 
many chemical species and a treatment of their interactions with a nuclear reaction 
network is necessary. We will leave such a treatment to future work (but see, e.g., 
H EH for discussions of such implementations). 



When using finite-temperature microphysical NSE EOS such as the LS EOS in 
GR hydrodynamics codes, two additional caveats need to be taken into account: 
(1) The thermodynamic potential from which all dependent variables arc derived 
is the Helmholtz free energy F, This makes the EOS a function of {p,T,Y e } while 
GR hydrodynamics codes such as GRID operate on the primitive thermodynamic and 
compositional variables {p, e, Y e }. Hence, in a typical EOS call it is first necessary to 
determine T{p, e, Y e ) through a root-finding procedure, before the dependent variables 
can be obtained through tri-linear interpolation in {p,T,Y e }. (2) In contrast to 
Newtonian hydrodynamics that involves only differences of the specific internal energy 
e, GR codes depend directly on e through its contribution to the matter stress-energy 
tensor. Hence, it is important to find and use a physically correct energy zero point 
and ensure that there are no rest-mass contributions included in e. 
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3.3. HShen EOS 



The HShen EOS [6l|, 162J is based on a relativistic mean-field model for nuclear 
interactions, assumes NSE, and is extended with the Thomas-Fermi approximation 
to describe the homogeneous phase of matter as well as the inhomogeneous matter 
composition. K of the HShen EOS is 281 MeV and the symmetry energy S v has a 
value of 36.9 MeV. The authors of the HShen EOS provide the baryonic componcnt[jj] 
in tabulated form only. The provided table is not uniformly spaced and has too low 
resolution to be used directly with fast tri-linear interpolation in simulations. Hence, 
we generate a finer uniformly-spaced table that has 18 points per decade in log 10 p 
from 10 3 - 10 15 - 36 gcm~ 3 , 41 points per decade in log 10 T from 10~ 2 - 10 24 MeV, 
and 50 points in Y e covering the interval 0.015 — 0.56. We interpolate all dependent 
variables from the original HShen table using the cubic Hermitc interpolation function 
given in [86[ modified to have monotonic interpolation behavior according to the 
prescription of [87| • The interpolation is performed first bicubic in p 1 T, then cubic in 
Y e . Alternatively to the just described, one could interpolate the Helmholtz free energy 
F and re-derive dependent variables by taking derivatives of F on the interpolated 
table (see, e.g., (8(|). We decided against this approach, since it would require quintic 
interpolation and the knowledge of the second derivatives of F at each point in the 
original table, some of which would have to be computed by taking second derivatives 
in the coarse original table. Also, compositional information cannot be obtained 
directly from F and would have to be interpolated from the original table. 

We perform the described interpolation at densities above 10 7 gcm -3 . For points 
with T > 100 MeV and T < 0.1 MeV we extrapolate most variables linearly, keeping 
only the compositions fixed. We add photons and electrons after interpolation using 
the routines of the Timmes EOS. At densities below 10 7 gcm -3 , we employ the 
Timmes EOS in the same fashion as described in the above for the LS EOS. 

We compute the maximum cold neutron star masses for the HShen EOS in the 
same way as for the LS EOS and find 2.24 M and 2.61 Mq, for the gravitational 
and baryonic value, respectively. The coordinate radius of the corresponding star is 
12.6 km. 



4. Neutrino Leakage and Heating 

4-1. Deleptonization and Electron Capture in the Collapse Phase 

Electron capture on free and bound protons leads to the emission of neutrinos that 
stream away from the core and carry away net lepton number at densities below 
~ 10 12 gcm -3 . Hence, one speaks of the deleptonization of the core. The change of 
the electron fraction Y e in the collapse phase due to deleptonization has important 
dynamical consequences. A reduction of Y e leads to a decrease of the mass of the 
homologously collapsing inner core whose kinetic energy is initially imparted on the 
supernova shock and which turns into the PNS core after bounce [lj . We take electron 
capture in collapse into account in GRID with the approximate scheme of Licbendorfcr 
[631 ] who observed that Y e of infalling mass elements depends primarily on the local 
matter density p and can be parameterized with rather high precision on the basis of 
radiation- hydrodynamic calculations. 
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Operator-split, after a hydrodynamics update, we compute the change in Y e , 

AY e = min[0,F e (p)-y e ] , (21) 

which ensures for consistency that a c hang e in Y e is either negative or 0. We use 
for Y e (p) the fitting formula given in [63J with parameters pi = 3 x 10 7 g cm~ 3 , 
p 2 = 2 x 10 13 gcm" 3 , Yi = 0.5, Y 2 = 0.278, and Y c = 0.035 corresponding to the 15- 
Mq model of 88|, evolved as model G15 by GRID also contains an interpolation 



routine to use numerical Y e (p) data. 

Electron capture leads to a change in the entropy (s, the specific entropy in units of 
fce/baryon) that is carried away by neutrinos leaving the core at densities below an 
assumed trapping density ptrap = 2 x 10 12 gcm~ 3 . The entropy change is given by 

As = -AY e ^-\ + ^- Ev . (22) 

E v is the energy of the escaping neutrinos (set to lOMeV). p p , p n , and p e are the 
proton, neutron, and electron chemical potentials including rest mass, respectively. 



Following [631 ] . we set As = if p p + p n + p e < E v and above ptrap- After updating 
the entropy, we use the EOS to update the specific internal energy e for consistency 
with the new Y e and s. 

We employ the outlined deleptonization scheme until core bounce (defined as the 
time when the peak entropy of the inner core surpasses 3 fes/baryon) and until 5ms 
after bounce for yet unshocked regions of the outer core that will settle in the high- 
density outer PNS and only in this way assume realistic postbounce Ye. 

4-2. Postbounce Deleptonization and Neutrino Heating/Cooling 

At core bounce a strong hydrodynamic shock wave is generated that travels outward 
into the outer core, heating and dissociating infalling heavy nuclei into nucleons. 
Electron capture occurs rapidly on free protons and a sea of electron neutrinos 
{v e ) builds up and is released in the v e burst when the shock breaks through the 
ncutrinospher deleptonizing the postshock region and leaving behind a "trough" 



in the Y e profile (e.g., |41(). The softening of the EOS due to dissociation of nuclei 
and postshock energy loss to escaping neutrinos lead the shock to stall and turn into 
an accretion shock soon after bounce. In the hot postshock region, electrons are less 
degenerate and positrons appear and are captured on neutrons, leading to a rise of 
the v e luminosity. In addition, in the PNS and in the postshock region, neutrinos and 
antineutrinos of all flavors are emitted by thermal processes. 

The simple Y e {p) parameterization discussed in the previous section 14.11 is 
not adequate to capture these effects and, in principle, a full neutrino energy- 
dependent radiation-hydrodynamics treatment would be needed for accurately 
capturing postbounce neutrino effects. Such a treatment may be added in future 
versions of GRID. In the present version of GRID, we approximate postbounce neutrino 
transport by a gray (energy-averaged) neutrino leakage scheme augmented with a 
simple prescription for neutrino heating in the postshock region. This approach 
captures the most important qualitative aspects of the postbounce evolution well 
and, as we demonstrate in section T6.21 is sufficiently quantitatively accurate to make 

% The neutrino-sphere is the effective "decoupling" surface of neutrinos where the optical depth r v of 
the supernova matter is 2/3. Its position depends strongly on neutrino energy. 
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reliable predictions of the time of BH formation and the maximum PNS mass in failing 
core-collapse supcrnovac. 

Our implementation in GRID combines elements of the neutrino leakage schemes of 
Ruffcrt ct al. [§4[ and of Rosswog & Liebendorfer [65|] . We consider three neutrino 
species, z/ e , D e , and v x . In the latter, we lump together fj, and r neutrinos and anti- 
neutrinos since they interact only by neutral-current processes in the core collapse 
context and have very similar cross sections. The mean (energy-averaged) optical 
depth is 



r„,(r) = / Kt(vi)Xdr , (23) 

J r 

where Kt(^i) is the mean transport opacity equal to the sum of absorptive and 
scattering opacitietQ for neutrino species v^. We follow [64[ in the calculation of K t (yi) 
and of the approximate neutrino degeneracy parameters [r\ Vi = fi^/T). We consider 
opacity contributions from neutrino scattering on neutrons, protons, and heavy nuclei 
and absorption of neutrinos (anti- neutrinos) on neutrons (protons). For heavy- lepton 
neutrinos that are never degenerate, we set r] Ux = 0. 77^ is known (1) in /3-equilibrium 
where ?/£ q = ?/ e + r) p — rj n (where we assume that the chemical potentials include rest 
mass terms) and (2) in the free streaming limit, where ^tream = 0. Furthermore, 
rfj^ = —r] 1 ^. In between the two regimes, the neutrino distribution function cannot be 
derived from first principles and neutrino transport is necessary for a correct estimate 
of 77„ e and -q^^. As an approximation, we interpolate between (1) and (2) using the 
optical depth, 

Th* =7£?(l-e- T "< ( ""< ) ) • (24) 

Note that r Vi depends on r\ Vi and vice versa. Hence, we iterate their calculation until 
convergence is reachecQ. 



Knowing r Vi and rj Vi , we use the leakage scheme of [651 ] to calculate the neutrino 
emission rates for the capture processes p + e~ — > v e + n and e + + n — > v e + p and 
thermal emission via electron-positron annihilation and plasmon decay to vv pairs. 



We modify the scheme of [65( in the following ways: (i) we use the interpolated r\ Vi 
from above instead of the equilibrium values suggested in [65j . (ii) we increase their 
diffusion time scale by a factor of 2 to obtain more reasonable neutrino luminosity 
predictions, and Jtii) for simplicity, we use the analytic thermal emissivities from 
[64l | . Following [65(, we then interpolate the effective volumetric energy loss Q^ k 
(erg/cm 3 /s) and effective number loss R^g (#/cm 3 /s) between the limits of diffusive 
emission (subscript "diff" ) and free emission (subscript "loc" ) using 

Xcff", = xttJi 1 + x!ot, /*X< ) . ( 25 ) 

where \ = Q for energy loss and \ = R for number loss (see [iiH ] for definitions and 
details). We define the neutrino luminosity seen by an observer at rest at radius r 
in the coordinate frame by summing up the effective energy emission rates from each 

+ Note that the opacities for neutrino number and neutrino energy transport differ. Hence, the 
optical depths for number and energy transport must be computed separately [64j . We neglect this 
subtlety and use the optical depths for energy transport throughout GRID. 

* Initially we choose K Vi (r) = 10 — 5 cm — 1 determine T Vi through l|23|l and iterate il 2 4 1) . For all 
subsequent times we use the previously determined value of t v . as a starting point, convergence 
(fractional difference in K Vi < 10~ 10 ) is typically reached after three iterations. 
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zone interior to r, transforming from the fluid rest frame (FRF) to the coordinate 
frame (CF), and applying the rcdshift (see |Appcndix B| for details), 





~a(r')~ 


f 

Jo 


a(r) 



(26) 



U^{r')[a{r')W(r'){\ + v{r'))]X(r')r' z dr' . 
For an observer at rest at r — oo (q(oo) = 1), 

L„ 4 (oo) =4tt / a(r') Q oS ^(r')[a{r')W(r'){l + v(r'))]X(r')r' 2 dr' . 
Jo 

It is useful to note the neutrino luminosity as seen by an observer at rest in the fluid 
rest frame at radius r, 



(27) 



a(r)W(r)(l+v(r)) ' 



(28) 



where the denominator transforms the luminosity from the frame of an observer at 
rest in the coordinate frame (l26l) to the fluid rest frame. 



Jf.,2.1. Neutrino Heating. In addition to the above leakage scheme, we include a 
parameterized heating scheme to mimick neutrino absorption in the postshock region. 
Heating occurs at intermediate to low optical depths where neutrinos begin to decouple 
from matter and a net energy transfer from neutrinos to the fluid is possible (see, e.g., 
(89j). The dominant heating processes are the charged-current capture reactions of v e 



on neutrons and v e on protons. We take the absorption cross sections from [65 



(l + 3g|) 

Theat,^ = "j <?Q-, 2\2\± ~ Je-f J \ Z ^J 

4 \TYl e C ) 

(l + 3 g j) (e 2 )g: 
OW.P. = j (m e c 2 ) 2 ^ ~ te+) ' ( ' 

where <tq is a reference weak-interaction cross section equal to 1.76 x 10 _44 cm 2 , 
qa ~ —1.25, and the Fermi blocking factors (1 — /j) are defined analogously to 
[64], Ir35| . In the postshock region the positron blocking term is negligible but the 
electron blocking term can be significant around the time of bounce. Following (89| , 
we set the mean squared neutrino energy to (e 2 )" t s = T(r Vi — |) 2 -7 r 5(?7" I s )/-7 : 3(?7" I s ), 
where T(r Vi = |) is the temperature at the neutrinosphere of species i, superscript ns 
denotes neutrinospheric values, and J-'niv) = Jo°° expfe^T+ i ^ s ^ e n * h Fermi integral 
(we approximate Fermi integrals via the formulae given in J90J ) ■ 

Given the neutrino luminosity L™ F (r) obtained from the leakage scheme (|28p. we 
write the local neutrino heating rate in units of ergcm _3 s _1 as 



(r) = /heat—/ o - CThcat,^ -X"i ( / e ^ > ( 31 ) 

where m u is atomic mass unit and the mass fraction = X„ in the case of v e 
absorption and Xj = X p for 9 e s. (l/F Ui ) is the mean inverse flux factor describing 
the degree of forward-peaking of the radiation field (e.g., 0, HH; (l/F Ui ) is 1 for 
free streaming and diverges at high optical depth). We estimate (l/F Vi ) by the 
interpolation (l/F Vi (r)) = 4.275r + 1.15, which reproduces the predicted values of 
4 at the neutrinosphere [89| and levels off at a value of 1.15 at low optical depth in 
the outer postshock region. We choose the latter value instead of 1, because (a) the 
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radiation field becomes fully forward peaked only outside the shock (e.g., [44J), and (6) 
the linear interpolation in r drops off too quickly compared to full simulations [44|, 
hence the higher floor value to compensate. Finally, we introduce the attenuation 
factor e~ 2T "i to cut off heating near and below the ncutrinosphere and the scaling 
factor /heat to allow for an ad-hoc increase of the heating rate. Once the heating 
rate for a computational cell is computed, we reduce the outgoing luminosity by the 
deposited power for overall energy conservation. In the coordinate frame (|26|) now 
becomes, 



f 


'a(r')' 


10 


a(r) 



[Qen^(r')-Q^(r')] [a(r')W(r')(l+v(r'))]X(r'y 2 dr' .(32) 



Along with the energy deposition goes a change in Y e which can be written as 

Qheat /Ohcat 

Ky - ~W"W (33) 

where we approximate the mean neutrino energies based on their neutrinospheric 
values as (e" s ) = T{r Ve = l)Mv™)/Mv™) M- 

To caution the reader, we point out that the simple gray heating scheme presented 
in the above is not self-consistent and cannot replace a radiation transport treatment 
that allows emission and absorption to balance. While we find that the combination of 
gray leakage/heating reproduces the overall qualitative dynamical features observed 
in postbounce radiation-hydrodynamic simulations, quantitative aspects are not 
captured as well. This is true in particular in highly dynamical situations shortly after 
bounce when we observe an unphysical rise of the electron fraction due to heating in 
the lower postshock region. 

We couple the neutrino leakage/heating scheme with the GR hydrodynamics in GRID 
through source/sink terms on the RHS of the GR hydrodynamics equations in MoL. 
Neutrino-matter interactions occur in the fluid rest frame where the total energy and 
number changes are given by 

r)0 /ohcat y^lcak oO ohcat i pleak {1A\ 

E — Vtotal — ^eff, total > U Y e ~ -"total + ^eff.total ) K' 3 *) 

where Q t h c t a a 'i and Qrfrjtotai are always positive or zero and Eg* and R^% tal may be 
positive or negative. Following [4jj, |56[, transforming these terms to the coordinate 
frame via the methods laid out in |Appcndix A[ we obtain the neutrino heating/cooling 
and deleptonization source/sink terms for the RHS in the MoL integration, 

= ctXRl , Q"^ = ctvWQl , Q^ E = oWQ% . (35) 



4-3. Neutrino Pressure 

Electron neutrinos above trapping density in the inner core during the final phases of 
collapse and in the postbounce PNS contribute to both the pressure and the specific 
energy density (with relative importance of up to ~ 10% around core bounce (9lt|). 
We neglect neutrino contributions to pressure and energy below ptrap where they are 
small, but otherwise follow [631 ] and assume electron neutrinos and antineutrinos to 
be a perfect Fermi gas. The pressure is then given by 

4-7T 
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where r\ v = pt, v /T and p v = fi e — p n + p p , where the chemical potentials include 
rest mass contributions. F 3 is the 3 rd Fermi integral which we approximate following 
92[. The specific internal energy of a relativistic Fermi gas of neutrinos is simply 
e v = ZP v /p. 

We treat neutrinos and fluid separately from each other and treat momentum transfer 
between the neutrino radiation field and the fluid approximately using the radial 
gradient of the neutrino pressure as suggested by (§3} . We couple this radiation stress 
into GRID's MoL integration of the GR momentum (S r ) and energy (r) equations via 
source terms (see |Appcndix A for a derivation; we neglect rotational effects in these 
source terms), 



-aW 



dP v 



Qr' M = -aWv^ 



(37) 



dr dr 

In addition to the force on the fluid due to the neutrino pressure gradient, we take 
into account the energy and "pressure" of the neutrino radiation field by adding P v 
and e v through the terms and r$ in d3| and Eqs. ([5]) and (jTS"]) . These contributions 
are derived by modifying the stress-energy tensor, 

'P + Pu ' 



T 



af) _ 



P 



1 + (e + e v ) + 



P 



g a P{P + P v ) 



and t£ are then given by [35[ 

T v m = P W 2 (e u + P v /p) - P v = {AW 2 - 1)P V , 
tI = pW 2 v 2 {e u + P v /p) +P„= {AW 2 v 2 + l)P v . 
We note that if rotation is included, v 2 in p0|) is replaced with v 2 



(38) 

(39) 
(40) 



lv 2 

3 "<P 



5. Code Tests 



In the following, we provide results from a set of standard and stringent relativistic 
hydrodynamics code tests for which analytic results exist. These involve two planar 
shocktube problems in section 15.11 the spherical Sedov blast wave problem in 
section I5T21 and Oppcnhcimer-Snyder collapse in section IS~3l Finally, in section [531 
we present results from a collapse simulation of a n = 3 polytrope and demonstrate 
convergence of the hydrodynamics scheme in GRID. With this selection, we test a broad 
range of aspects of potential problems to be addressed with GRID: special relativistic 
effects, geometrical effects, and fully general-relativistic collapse dynamics. 



5.1. Relativistic Shocktube 

We assume flat space and planar geometry and perform the two relativistic shocktube 



tests proposed by [93j. We use a T— law EOS with T = 5/3 and a grid of length 1 
with a cell spacing of dx = 0.001. The starting values of the density, pressure and 
velocity are summarized in Table [TJ The left panel of figure [1] shows the exact results 
for velocity, density, and pressure of the mildly-relativistic problem #1 at t = 0.4. 
Superposed are the numerical results obtained with GRID that reproduce the exact 
results nearly perfectly. Problem ^2 is a more stringent test and involves Lorentz 
factors of up to 6 in the forward propagating shock and a very thin shell of trailing 
matter. As shown in the right panel of section [TJ GRID reproduces the exact solution 
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Table 1. Initial conditions for two relativistic shocktube problems as presented 



at t = 0.4 very well almost everywhere, but fails to completely resolve the thin shell 
of relativistic matter. This is most likely due to the rather diffusive nature of the 
HLLE Riemann solver employed in GRID (see, e.g., [9l|, [94[ for comparable results 
obtained with a nominally more accurate scheme). In an attempt to obtain results 
closer to the analytic solution we use 3 rd order Runge-Kutta time integration for this 
test case. These deviations are not worrying since the shocks that obtain in stellar 
collapse are much less relativistic than that of problem #2. If GRID were to be applied 
to ultrarelativistic outflows (e.g., in a GRB), a more precise treatment of the Riemann 
problem would likely be necessary. 

5.2. Sedov Blast Wave 

The above shocktube tests demonstrated the ability of GRID to capture shocks and 
solve the special-relativistic hydrodynamic equations in planar geometry. Here we 
go back to Newtonian hydrodynamics and test instead spherical hydrodynamics 
with Sedov's blast wave problem [95| . For a comparison with a large number 
of hydrodynamics codes, we use the initial conditions of [96j |. The grid setup is 



Problem #1 
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Figure 1. Relativistic shocktube simulations: Initial conditions taken from [93j j 
and provided in Table[T] The pressure, density, and velocity are shown at t = 0.4 
for problem #1 (left panel) and problem #2 (right panel). For reference, in 
both figures the pressure is denoted by boxes (red), density by circles (blue) and 
velocity by diamonds (green). The analytic solution is denoted by the solid line. 
Both problems were run with a Courant factor of 0.5 and 3 rd order Runge-Kutta 
integration. 
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Radius 



Figure 2. The Sedov blast wave problem and exact solution at t = 0.1. Shown 
are the numerical results with the exact solution underlying the various curves of 
density, pressure and velocity. Both the exact solution and the numerical result 
are normalized to the analytic value at the shock. p s = 4, P s = 252.255 and 
v s = 13.757. 



in spherical geometry with (dimcnsionless) r max = 10 and N = 400 cells which 
corresponds to the maximum mesh refinement level used in (9(| . We deposit a constant 
specific internal energy into a sphere of radius r = 0.0875, corresponding to a total 
(dimcnsionless) energy of E a = 10 5 , into a background medium of (dimensionless) 
Po = 1. We set the background energy density to an insignificant amount and use a 
T-law EOS with T = 5/3. Figure [2] depicts the comparison of our numerical solution 
with the exact result for density, velocity and pressure at t = 0.1 normalized in such 
a way that the value of all variables at the shock is 1. GRID performs very well in the 
region behind the shock and provides an adequate, though not perfect, solution near 
the shock. 

In addition to the Newtonian Sedov blast wave problem, we have also considered its 
relativistic variant discussed in (9?]]. These authors used 17 levels of adaptive mesh 
refinement (AMR) and we find that the lack of AMR in GRID makes it computationally 
impossible to adequately resolve the relativistic Sedov problem. This, however, is not 
a problem for the application of GRID to the stellar collapse problem, since the shocks 
appearing there are only mildly relativistic. 



5.3. Oppenheimer- Snyder Collapse 

For the final test problem for which an exact solution exists, we perform a simulation 
of the Oppenheimer-Snyder collapse (OSC) [98| of a constant-density sphere of 
pressureless (P = 0) dust. The exact solution of OSC in RGPS spacetime has been 
laid out by jo! [l00jj. We choose M = M , = 1OM . We perform the OSC test 
with the standard version of GRID described in section[2]of this paper and do not make 
special adjustments for the code to operate with P = 0. Hence, we set the pressure 
to a small, but non-zero value, using a polytropic EOS with K = 10 -20 and T = 5/3. 
In the artificial atmosphere outside the dust ball, we set the density to lgcm -3 . We 
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Figure 3. Oppcnhcimcr-Snydcr collapse of a pressureless dust ball. Shown are 
the numerical (plus symbols) and exact (solid lines) density (left panel) and lapse 
(right panel) profiles for various times. The density is normalized to the density 
at t = 0. The simulation uses 9000 equally spaced grid points across the domain 
of 20 Mq. Initially one solar mass is distributed with constant density in a sphere 
of radius 10 Mq. For clarity, we show only every third data point. 



use 9000 equidistant zones to model OSC with GRID. 

In figure |3l we compare numerical and exact density and lapse profiles of OSC at 
t = 30, 35, 40, 43 and 60 Mq. Following [33 |. we normalize the central density to 
the value at t = 0. The overall agreement is excellent. However, we notice two slight 
deviations: (1), near the origin, we observe a small build up of material. This is 



present also in the OSC test of [34| and probably due to diverging terms near the 
origin. We do not notice this effect in our stellar collapse calculations, most likely 
because of the stabilizing effect of the large pressure in the PNS. (2), at late times 
(t > 50Mq), the numerical a decreases more slowly then its exact counterpart and 
begins to deviate significantly at a(r = 0) < 0.001. We attribute this to numerical 
inaccuracies developing due (a) to the metric coefficient X becoming singular as 
Ri, —> 2M®, (6) to the extreme density gradient developing at the surface at late 
times, and (c) to the fact that we use the standard version of GRID without special 
adjustments for the OSC problem (as, e.g., made by [IH). 

5.4- Hybrid Core Collapse: Convergence 

In this section, we present simulations of nonrotating core collapse and present proof 
of convergence for GRID. We utilize the hybrid EOS described in section |3"7T1 taking 
P*! = 1.28, T 2 = 2.5, r th = 1.5 and K = 4.935 x 10 14 [cgs]. Following j9l|, we use as 
initial data an n = 3 polytrope with a central density of p c = 5 x 10 10 g/cm 3 and a K 
value as above and initially zero radial velocity. We simulate the evolution with GRID 
for equally spaced grids of three different resolutions (-/V ZO ncs = 500, 1500 and 4500) 
to test the self- convergence of the code. The self- convergence factor at convergence 
order n of a quantity q is given by, 

n gi - (12 (tfaiT - idx 2 ) n 

V 92-93 {dx 2 ) n - (dx 3 ) n ' 1 1 
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where is the numerical result from the simulation with the corresponding resolution 
and dx is the zone width. For this convergence test, dxi = 3dx 2 — 9dx^. In the 
lower panel of figure HI we show the self-convergence of M grav at t = —3.3 ms (before 
bounce) as well as at t — 16.6 ms and t — 26.6 ms after bounce. 




Radius [km] 



Figure 4. Radial density profiles and self-convergence for core collapse using 
the hybrid EOS. Top: Density profiles of the core collapse for various times 
including in the prcbouncc phase, and after the shock has propagated through 
~ 300 and 600 km. We show the low resolution profile (segmented lines) as 
well as the high resolution profile (solid lines) for comparison. Bottom: Self- 
convergence of the enclosed gravitational mass, m(r). Dotted lines at Q=3 & 9 
denote expected values for 1 st and 2 nd order convergence. 



We generally see the expected 2 nd order convergence (Q=9) in smooth parts of 
the flow, but note several interesting features: (1) before bounce (red, dot-dashed 
curve) and near 120 km where the convergence spikes, the velocity is peaking, causing 
a reduction in convergence. (2), during the postbounce phase, convergence in the 
shocked region drops to 1 st order, this is characteristic of HRSC schemes in the 
presence of shocks. (5), finally, during the postbounce phase for r < 20 km, the 
steepness of the density gradient at the PNS surface and the coarseness of the grid lead 
to local non-convergence. We note that the lowest resolution used here is dx ~ 2 km 
and that deviations in the density profile compared to higher-resolution simulations 
can be seen in the top panel of figure 2] 



6. Sample Results for a 4O-M Star 



In the following simulations we use the single-star, non-rotating Mzams = 40 M, 



0! 



solar-metallicity presupernova model of Woosley & Weaver [88( (model s40WW95 
hereafter). This model has an iron core mass of 1.98 Mq. We set up a grid of 1000 
zones that is logarithmically spaced from r = 20 km outward, extending to a radius of 
1.15 x 10 5 km where the density drops to 200gcm~ 3 . There is 14.7 Mq of baryonic 
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Figure 5. Black hole formation with rotation and the hybrid EOS. Left panel: 
Central densities for various initial angular velocities (top) and fractional error of 
the conserved quantities M grav and J (bottom). Right panel: T/\W gr av | near 
bounce and (inset) over the entire simulation. f2(r) is set through {42}. 



material within this density cutoff. Inside r = 20 km, we use an equidistant grid with 
a spacing of 100 m. Such high resolution is necessary to resolve steep gradients at 
the PNS surface at late times (t > 0.5 s). Near the origin, we increase the zone size 
gradually to ~ 700 m for improved stability but for rotating runs we find it necessary 
to maintain the fine grid spacing all the way to the origin to capture the correct 
angular velocity profile. 



6.1. Rotating Core Collapse and Black Hole Formation in a 40-Mq Star using the 
Hybrid EOS 

To show the effects of including rotation and to further demonstrate the use and 
usefulness of the hybrid EOS (see section I3.1[) for exploratory studies, we perform 
a set of collapse simulations to black hole formation. We set Y\ = 1.30, T2 = 2.5, 



r t h = 1-34 and impose rotation according to the rotation law (see, e.g., [771. l80j ) 



• - 1 

n(r) = e— 
w no 







1 + 









rads" 1 , (42) 



where we vary £ from to 5 and A is a parameter governing the degree of differential 
rotation. We choose A = 1000 km which leads to roughly uniform rotation within 
the inner core as predicted by stellar evolutionary calculations (e.g., [zij]). As an 
additional test of GRID, we show in the lower part of the left panel of figure [5] the 
relative error in total angular momentum and gravitational mass in the most rapidly 
spinning simulation. GRID conserves angular momentum to better then one part in 
10 4 and M grav to one part in 10 6 until the onset of BH formation when the resolution 
becomes insufficient to fully resolve the huge gradients in the collapsing PNS. 

We show in the top part of figure [5] the evolution of the central density in the 
simulated models. Due to the choice of Ti, rotation has little influence on the 
prebounce dynamics [H^]. The hybrid EOS qualitatively captures the stiffening of 
the EOS at nuclear density that leads to core bounce. Owing to the small value of 
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r t h, the shock stalls soon after bounce and accretion on the PNS continues. Slowly 
spinning models accrete rapidly and collapse to a BH after only 200 ms. Centrifugal 
support becomes dynamically relevant in more rapidly spinning cases, decreasing 
the accretion rate and delaying BH formation. The right panel of figure [5] depicts 
the evolution of the rotation parameter T/|W grav |. Its systematic s are very similar 



to what has been observed in multi-D simulations (e.g., [53J, |77j, llOlj ). T/|W gra v| 
reaches a local maximum at bounce, then decreases as the PNS reaches its postbouncc 
quasi-equilibrium. New and not shown before is the evolution of T/|W g rav| near 
to BH formation. T/|Wg rav | increases only slowly after bounce (note that, in a 
calculation with neutrino transport or leakage, the postbounce T/|Wg rav | would 
increase faster jz3|), but near BH formation grows nearly exponentially during PNS 
collapse. Rotation, in particular when it is st rongly differential, can increase the 
maximum mass of the accreting PNS (e.g., [l02| ) . We finely BH birth masses of 1.89- 
1.97 Mq for the set of rotating hybrid-EOS models considered here. This increase in 
the maximum mass is modest, primarily because our PNS cores are rather uniformly 
spinning (in agreement with 53l. l77|). We point out that our present treatment does 
not consider angular momentum redistribution by multi-dimensional effects or effective 



viscosity which may be present in realistic systems (see, e.g., |54j, |76j and references 
therein) . 

Finally, we note that for the nonrotating (£ = 0), the evolution with GRID continues 
until a central value of the lapse function of 3 x 10~ 10 and a maximum value of 
y/g r r = X of ~ 21.1. These ar e ex cellent values in comparison to previous studies 
on BH formation in RGPS EJ llOOj ]. In the rotating case, the evolution terminates 



somewhat earlier due primarily to numerical issues near the origin at very large 



6.2. Nonrotating Collapse and Black Hole Formation with Neutrino Leakage/ Heating 
in a 40 -Mq Star 

In this section we show example results employing GRID's leakage/heating scheme and 
finite-temperature EOS. We use the s40WW95 progenitor and the LS180 EOSlB], Y e {p) 
parameterization prc-bounce, our standard leakage/heating scheme after bounce, and 
no rotation. We show results for both /heat = (losses only) and /heat = 1- In 
figure |6l we compare the shock radii of these two runs and neutrino luminosities of 
the /heat = 1 run (left panel) as well as the Y e radial profiles at 50 ms after bounce 
(right panel). We note that the total luminosity is L Ue + Lp s + AL V and is corrected 
for redshift through (|32[) with r = oo, but, nevertheless, is somewhat higher (up to 
~ 20%) than predict ed by full Boltzmann radiation-hydrodynamics calculations using 
the same progenitor |l03l . Il04j . The time until BH formation in the case of /heat = 1 is 
^bh =511 ms and the baryonic mass inside the shock of the last stable configuration is 
2.25 Mq. We compare this to two other studies of BH formation in ID with the same 
progenitor model and EOS, but with two different im plem entations of GR Boltzmann 
neutrino transport. These studies a re Fi scher et al. [103| who found £bh = 435.5 ms 
and 2.196 M© and Sumiyoshi et al. |l04j . who found £bh = 560 ms and 2.1 Mq. Our 
result is very close to these more accurate studies which gives us confidence in the 

J In RGPS, a coordinate singularity develops at R = 2M upon BH formation. We define here the BH 
mass to be M grav inside the radius that corresponds to the maximum X. This is an approximation 
and is subject to errors due to our finite resolution grid. 

ft The lower bound on our EOS tables is 1000 gem -3 , we bring the outer boundary into p = 
2000 g cm -3 for this example. 




Figure 6. Left panel: Shock radius (thick line, right ordinate) and neutrino 
luminosities (thin lines, left ordinate) as a function of postbounce time in a 
nonrotating leakage+heating (/heat = 1) simulation with the 4OM0-model of [88| 
run with the LS180 EOS. Shown also is the shock radius evolution (dashed thick 
lines) in a simulation without heating /h ca t = 0- Right panel: Y e profiles of both 
simulations at 50 ms after bounce, corresponding to the maximum shock radius 
of the fheat = 1 simulation. Shock radii, electron and anti-electron neutrino 
neutrinospheres are marked for both the fheat = 1 an d /heat = simulations. 



robustness of the heating/leakage scheme in GRID. 

The right panel of figure [6] depicts the Y e profiles at 50 ms after bounce. The 
characteristic trough in Y e behind the shock is captured by our leakage/heating 
scheme, but we find that our simple heating scheme converts too many of the postshock 
neutrons back to protons at early times, leading to too high values of Y e in the lower 
postshock region between ~ 30 — 60 km. 

To conclude this section, we note that, due to the computational efficiency of our 
scheme, each of our simulations took only ~ 6 CPU hours from iron core collapse 
through BH formation on one core of an Intel Xeon X5550 (Nchalem) machine. 



7. Summary and Concluding Remarks 



In this paper, we have presented the details of our new open-source Eulerian 1.5D GR 
hydrodynamics code GRID. GRID is intended primarily for the simulation of stellar 
collapse to neutron stars and black holes and, for the first time in the ID GR 
context, includes an approximate way of accounting for stellar rotation consistent 
with that used in state-of-the-art calculations of stellar evolution (e.g., 74|). Using 
this scheme, we have presented rotating long-term postbounce simulations towards 
black hole formation using a 40- M Q supernova progenitor model and showed how the 
simple analytic hybrid EOS can be used to capture many qualitative aspects of this 
phenomenon. 

As we have demonstrated in this paper, GRID performs well in standard tests and, 
despite its simplified neutrino leakage/heating scheme, still yields overall results in 
the case of failing core-collapse supcrnovae and black hole formation that measure 
up qualitatively and to some extent also quantitatively to those obtained with full 
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Boltzmann neutrino transport in ID Lagrangian codes 103l . 1104 ( . 

Many ID GR (radiation)-hydrodynamics formulations have been presented in the 
past ~50 years. Yet, there is presently no open-source ID GR stellar collapse code 
available to the broader community. The primary motivation driving the development 
of GRID is the need for such an open-source code that may be used as a codebase, 
benchmark, and testbed for improved modeling technology to be included in multi-D 
GR codes addressing core-collapse supernova explosions, but also failing core-collapse 
supernovae, black hole formation, and the post-merger evolution of binary neutron-star 
and neutron-star - black hole coalescence. Equipped with an approximate neutrino- 
leakage scheme to capture the key effects associated with neutrino heating and cooling, 
the version of GRID discussed in this paper is a solid starting point for the next 
generation of astrophysically-relevant multi-D GR simulations. 

The current limitations of GRID due to its gray leakage and simplified heating 
scheme are obvious. We will continue to develop and improve GRID and intend to 
include as a next step energy-dependent radiation transport in the multi-group flux- 
limited diffusion appro xima tion (MGFLD) and/or in the isotropic diffusion source 
approximation (IDS A, }l05l |). 
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Appendix A. Evolution Equation Derivation 

In this appendix we derive the evolution equations for the conserved variables 
D, DY e , S r , S$ and r used in GRID and presented in section |2~21 and |2~51 GRID uses the 
spherically symmetric metric g uv = diag(— a 2 , X 2 , r 2 , r 2 sin 2 8) with a = exp ($(r, t)) 
with <3?(r, t) defined through ([5]), X = (1 — 2m j r ' t - ) )~ 1 / 2 where m(r,t) is the enclosed 
gravitational mass at coordinate radius r. We assume the matter to be a perfect fluid 
described by a mass current density of J M = pu^ and a stress-energy tensor, = 
phu^u v + g^ u P where p is the rest mass density, P is the fluid pressure, h = 1 + e + P/p 
is the specific enthalpy with e the specific internal energy; u M = (W/a, Wv/X, 0, 0) is 
the fluid 4- velocity (without taking into account rotation) with W = 1/ \J\ — v 2 is the 
Lorentz factor and v is the physical radial velocity. 

While evaluating the covariant derivative of the stress-energy tensor and matter 
current density, we make use of the following formula, 

= -L (V=ffJ") u (A.l) 



and 



V^T"" = {V^gT^) ^ + r^T" Q , (A.2) 



-9 



where y/—g = aXr 2 is the determinant of the metric and Y u ' afi are Christoffel symbols 



a 1 1 

and arc defined through derivatives of the metric, 

T^QM = 7j9 Vl: 5 (5w9,a + 9aP,n ~ 9an,p) ■ ( A -3) 

For our metric, all non-zero Christoffels are given in Table [All r" au is symmetric in 
the last two indices, duplicates are omitted. 

It is useful to note the following derivatives needed in the derivation of the evolution 
equations: 

a r $(r, t) = X 2 ™ + Aixr(P + phW 2 v 2 ) , (A.4) 
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r* M = 




r r 6)6» — 


TF 1 


1 lr - 


d r (/>(r,t) 


i 00 = 


r sin 2 # 


r* - 


a- 2 ^d t m(r,t) 




1 

r 


1 tt — 




1 00 - 


— sin cos 9 


i tr - 


^-d t m(r,t) 




l 

r 




£(^ m ( r ,t)-2fcil) 


r% - 


cos 8 
sin 



Table Al. Connection coefficients. 



d r X 

d t X 
d r m 
d t m 



X 3 
X 



d r m m 



3 dtTn 
r 



Anr^phW 2 - P) , 
,aphW 2 v 



= - 47rr z 



X 



Appendix A.l. Source Terms 



and V M T^ 



(A.5) 

(A.6) 
(A.7) 

(A.8) 



0. Since we treat 



The evolution equations follow from V ^J^ 
neutrinos through a leakage scheme, we add in neutrino source terms explicitly to the 
RHS of these equations. The neutrino physics of GRID occurs in the rest frame of the 
fluid; in this frame the energy and lepton rates are calculated with the neutrino leakage 
scheme, Q° E and Ry are given in (f34|). Momentum exchange in the fluid rest frame 
is taken into account approximately via Q° AI = — ^f- where the gradient is evaluated 
numerically in the coordinate frame. This introduces a slight inconsistency, since in 
a full radiation-transport treatment the momentum transfer is computed fully locally 
via the second angular moment of the local neutrino radiation intensity [58| . 

By writing the evolution equations in the comoving orthonormal frame of the fluid 
(fluid rest frame, [FRF]) with 4-velocity u = (1,0,0, 0)frf and unit radial normal 
ft = (0, 1, 0, 0)frf and expressing them as frame-independent tensor equations we can 
derive expressions for the evolution equations in any frame. For the lepton fraction, 



d t {pY e ) = R° Yc , 

dtipY^) = Ry e i 
d»{pY K vP) =R Ye , 

V„(pYeu") = R Yc ■ (A.9) 

We write the energy and momentum source terms in the fluid rest frame as a 4- 
vector, q = (Q E , Q%' 0> 0)frf or m frame-independent notation, Q E u + Qj^n. In the 
fluid rest frame, the evolution equations for energy and momentum become, 

d t T u = Q% = q t , (A.10) 

and 

d t T tr = Q° M = q r , (A.ll) 
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or in frame-independent tensor notation, 

Vj.T"" = q" . (A.12) 

For the evolution equations, we must transform q from the fluid rest frame, 
to the coordinate frame (CF) of GRID. In a general frame n is a vector that 
is both i) normalized and it) orthogonal to u. In the CF of GRID, where u 
is the 4-velocity, these two conditions (along with the assumption of spherical 
symmetry) on ri give n = (Wv/a, W/X, 0, 0) CF . q in the CF then becomes q = 
(^■{Q% + vQ° M ), jf(vQ° E + Q%), 0, 0) cp . This can also be derived via a Lorentz 
transformation. In principle, non-zero rotation will give rise to source terms for 
the ^-momentum evolution through and modify the radial source terms q r . In 
consideration of the significant approximations already present in both our neutrino 
leakage scheme and in our treatment of rotation, we neglect the influence of rotation 
on the source terms. This is justified as long as v v <C c. 



Appendix A. 2. GRID Evolution Equations 



In the coordinate frame of GRID where u' 1 = (W/a, Wv/X, 0, 0), the continuity 
equation, V M J M = gives the evolution of the rest mass density, 

Vj.fow") =0, 



1 



d t (D) 



'-9 

a I 



d r 



-pWv 



X 



-Dv 



0. 



0. 



(A.13) 



The evolution of the electron fraction Y e follows a similar derivation but contains a 
source term from the neutrino leakage scheme. In the coordinate frame of GRID (IA.9I) 
becomes, 



1 



1 

~aX 



-pWY e 



I 



- P WY e v 



dt (XpWY e ) + —d r —XpWY e V 



1 



dt(DY e ) + —d r [ -r^-DY e v 



,.2 



X 



aXR Y . (A.14) 



The momentum evolution equation for GRID is obtained by evaluating (|A.12[) with 
v = r. 



V M T" r 



= q 



-9 T» r ) 



-9 q 



-9 , 



d t (phW 2 v) + -^d r ( — (phW 2 v 2 + P) ) = aXq r - aX(T r vt T tv + T r vr T rV 



+ 3$ ( ^ ( ' 



aXq r - aX(T r tt T tt + T r rt T tr 
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d t (S r 



ar 

~x 



{S r v + P) 



T r tr T rt + Y r rr T rr + TfoT** + T r 99 T ee ) 

= - aX 2 T^dtm 

r aX 

m\ ( phW 2 v 2 + P 



X^_ 

r 



2P a 2 phW 2 - P 



X 2 r X 2 



aXq r 



d t (S r ) + —d r 



OLT 

-^(S r v + P) 



aX 



{S r v-T-D) [8irrP + — 



Pm ^ 2P 



X 2 r 



+ aW{vQ° E + Q° M ). (A.15) 



where in the last step we have reorganized the source terms to the form of [34j using 
the derivatives defined in (|A.4[) - (|A.8|) . If non-zero, = Wv v /r leads to an additional 
term {aphW 2 v 2 v sin(9) 2 /Xr) arising through T^T^ on the RHS of EH averaging 
this term over the spherical shell gives 2/3 aphW 2 v 2 / Xr. When rotation is included, 
the evolution equation for = phW 2 v v r is, 



= 0, 



d t (aXr 2 g H T^) + d r (aXr 2 g^T^) 



9 r^T^ , 

— 9 {T^Tf+T*T;) , 



d t (XphW 2 Vip r) + \d r ( ^—phW 2 v v rvX 1=0, 



X 



phW 2 v v r 
X 



-d t X - —d r X 



X 



aphW 2 v ip vX [4irr 2 P 



(A.16) 



The energy evolution equation for GRID is derived by taking v = t in (|A.12[) . 



d t (^{phW 2 - P)\ + ±d r (^-phW 2 v^j = aXqt ~ aX (i^T* + I^T^ 



X 
a 



d t {T + D) + -d r 



1 „ ( ar 2 



X 



S r 



aXq* - aX{T\ t T u + 2Y\ r T 

X 

s a 

aphW 2 v /X 



T l rr T rr ) - (PhW 2 - P)d t 



X 



-Or 



REFERENCES 



30 



dt(r + D) + \d r ( ^-S r ) = a 2 q t 



1 „ / or 2 



9t(r) + 3 9 r l^-(S r -Dv)\ =aW(Q° E + vQ° M ). (A.17) 



where in the last step we use the continuity equation (|A.13[) to subtract out the 
evolution of the rest mass density, obtaining the evolution equation for r. A non-zero 
does not contribute source terms to this evolution equation. 

Appendix B. Neutrino Luminosities 

The luminosity computed from the neutrino leakage scheme is derived in the rest 
frame of the fluid. We require knowledge of the neutrino luminosity as measured by 
an observer at rest in the coordinate frame to determine i) the luminosity measured 
by an observer at rest at infinity and ii) the luminosity in the fluid rest frame at some 
other coordinate radius for our neutrino heating scheme. We derive these relationships 
by assuming the neutrinos are emitted radially in the fluid rest frame with energy 

£>FRF 

In the fluid rest frame (FRF) , the 4- momentum of the (massless) neutrino is p a = 
(£ frf ,£ frf ,0,0)frf- We use the orthonormal tetrad in | Appendix A.l in the fluid 



frame, u = cq = (1,0,0, 0)frf and n = e\ = (0, 1,0, 0)frf, in the coordinate frame 
(CF), vP = &q = {W/a,Wv/X, 0,0) CF and n 13 = e{ = (Wv/a, W/X, 0, 0) CF . In this 
we have neglected rotational effects which will be small for v v <C c. Transforming p a 
to the coordinate basis of GRID, 

\a X J CF 

An observer at rest in the coordinate frame (U a = (1,0, 0, 0)cf) then sees the neutrino 
with energy, 

E CF = -p- U - -g a a P p U a = a 2 E FRF — (l + v) = aW(l + v)E FRF . (B.2) 

a 



Noting that (see jl06l |. eq. 25.25), for massless particles emitted from rest at r and 
observed by a observer at rest at r', A(r)|(7 o( r )| -1 ^ 2 = M r ')ISoo( r ')l -1 ^ 2 implies, 

£°V) _ \_ _ IffooWl 1 / 2 _ a{r) 

E CF (r) \ r - |.goo(r')| 1/2 a(r>) ' 1 ^ 

this is the redshift formula for particles leaving a gravitational well. 



